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ABSTRACT 

The mechanisms that maintain thermal balance in the intracluster medium (ICM) 
and produce the observed spatial distribution of the plasma density and tempera- 
ture in galaxy clusters remain a subject of debate. We present results from numerical 
simulations of the cooling-core cluster A2199 produced by the two-dimensional (2-D) 
resistive magnetohydrodynamics (MHD) code MACH2. In our simulations we explore 
the effect of anisotropic thermal conduction on the energy balance of the system. The 
results from idealized cases in 2-D axisymmetric geometry underscore the importance 
of the initial plasma density in ICM simulations, especially the near-core values since 
the radiation cooling rate is proportional to n e 2 . Heat conduction is found to be non- 
effective in preventing catastrophic cooling in this cluster. In addition we performed 
2-D planar MHD simulations starting from initial conditions deliberately violating 
both thermal balance and hydrostatic equilibrium in the ICM, to assess contributions 
of the convective terms in the energy balance of the system against anisotropic thermal 
conduction. We find that in this case work done by the pressure on the plasma can 
dominate the early evolution of the internal energy over anisotropic thermal conduc- 
tion in the presence of subsonic flows, thereby reducing the impact of the magnetic 
field. Deviations from hydrostatic equilibrium near the cluster core may be associated 
with transient activity of a central active galactic nucleus and/or remnant dynamical 
activity in the ICM and warrant further study in three dimensions. 

Key words: galaxy clusters, intracluster medium, magnetohydrodynamic simula- 
tions. 



1 INTRODUCTION strongly suggesting that one or more heating mechanisms 

must be balancing radiative cooling in galaxy clusters. 

Clusters or galaxies are the largest structures m the uni- 
verse. The majority of their barionic mass is in the form of A wide variety of heating mechanisms have been inves- 
a hot plasma filling the space between galaxies, the intra- tigated, including i) heat transport from the outer regions 
cluster medium (ICM). The cooling time of the ICM near of the cluster to the central cooling gas by conduction (e.g. 
the center of the clusters due to radiative losses via X-ray Binney & Cowie 1981; Tucker & Rosner 1983; Bertschinger 
emission is sh o rter than the Hu bble time for most clusters & Meiksin 1986; Bregman & David 1988; Gaetz 1989; Rosner 
l|Sarazinlll98"6l ; lEdge et al.lll992l ; Peres et al. 1998; Sander- & Tucker 1989; David et al. 1992; Pistinner & Shaviv 1996; 
son et al. 2006). In the absence of heating the short cooling Narayan & Medvedev 2001; Voigt et al. 2002; Fabian et al. 
times (as short as ~0.1 Gyr) are expected to lead to cooling 2002; Zakamska & Narayan 2003; Kim & Narayan 2003a; 
of the ICM plasma and its inflow towards the center of the Markevitch et al. 2003; Dolag et al. 2004; Xiang et al. 2007; 
cluster at a rate that was estimated to be as large as ~ 10 3 Pope et al. 2006; Balbus & Reynolds 2008; Bogdanovic et 
M yr' 1 (Cowie & Binney 1977; Fabian & Nulsen 1977; al. 2009; Parrish et al. 2009; Parrish et al. 2010); ii) energy 
Mathews & Bregman 1978; for a review, see Fabian 1994). injection from active galactic nuclei (AGN) via jets, bubbles, 
However, high resolution X-ray spectroscopy from XMM- radiation and sound waves (e.g. Binney & Tabor 1995; Ciotti 
Newton and Chandra have failed to find evidence of such & Ostriker 2001; Briiggen & Kaiser 2002; Omma et al. 2004; 
"cooling flows" (for a review see Peterson & Fabian 2006), Ruszkowski et al. 2004a,b; Voit & Donahue 2005; Fabian et 

al. 2005; Nusser et al. 2006; Mathews et al. 2006; Chandran 
& Rasera 2007; McCarthy et al. 2008; Briiggen & Scanna- 

* E-mail: Ioannis.G.Mikellides@jpl.nasa.gov pieco 2009); iii) heating due to galaxy motions via dynami- 



2 /. G. Mikellides, K. Tassis and H. W. Yorke 



cal friction (e.g. Schipper 1974; Miller 1986; Just et al. 1990; 
Fabian 2003; El-Zant et al. 2004; Kim et al. 2005; Dekel & 
Birnboim 2008); iv) turbulent mixing (e.g. Loewenstein & 
Fabian 1990; Deiss & Just 1996; Ricker & Sarazin 2001; Kim 
& Narayan 2003b; Cho et al. 2003; Dennis & Chandran 2005; 
ZuHone et al. 2009); v) magnetic field reconnection (Soker 
& Sarazin 1990) ; vi) viscous heat ing regulated by pressure 
anisotropy (e.g. iKunz et al.ll2010l ); or vii) a combination of 
these (e.g. Ruszkowski & Begelman 2002; Brighenti & Math- 
ews 2002; Brighenti & Mathews 2003; Voigt & Fabian 2004; 
Fujita & Suzuki 2005; Pope et al. 2006; Conroy & Ostriker 
2008; Guo et al. 2008). 

Due to the extremely high Hall parameter of ICM plas- 
mas the impact of the magnetic field on the thermal bal- 
ance of the ICM has been of particular interest; exten- 
sive multi-dimensional magnetohydrodynamic (MHD) sim- 
ulations have been emplo yed in the last seve r al ye ars to 
assess anisotropic effects. IParrish fc Quataertl (|2008) per- 
formed two-dimension al (2-D) and 3-D simulations with 
the M HD code Athena l|Gardiner fc Stonell2008l ; IStone et~ai1 
|200S| ) focusing on the evolution and sat uration of the h eat- 
flux-driven buoyancy instability ( HBI) (jQua taert 2008) for 
various magnetic field strengths. IBogdanovic et aL I (|2009l ) 
also performed 3-D simulations with Athena. They carried 
out pre-simulations, whereby starting with an isothermal 
ICM they allowed radiative cooling without thermal con- 
duction to obtain a cooling core state, and then used it as 
the "initial condition" for their self-consistent MHD simula- 
tions. For the initial magnetic field the authors' numerical 
experiments included a case of rotational field lines. They 
concluded that thermal conduction alone cannot be the so- 
lution to the cooling flow problem. Moreover, they hypoth- 
esized that the nature of the AGN feedback is one of "stir- 
ring" rather than heating, that is, AGN activity periodically 
disrupting the azimuthal field structure allowing thermal 
conduction to sporadically heat the core. The level of ef- 
fectiveness of "stirring" however remains an open question. 
For example, cooling-core clusters like the A2199 contain 
prominent metal gradients that could preven t violent pro- 
cesses from stirring the ICM. In a similar studv lParrish et al.l 
(2009) performed a wide range of numerical studies of the 
A2199 cluster with tangled and radial initial magnetic field 
topologies and str engths of 1 nG to 1/xG, extending the 
parameter study of IBogdanovic et alj (|2009h . They demon- 
strated that the effective reduction of the anisotropic heat 
flux was less than 10% of the Spitzer value due to a re- 
orientation of the magnetic field lines relative to the radial 
direction induced by the HBI. This reduction was found to 
occur withing several times the HBI growth time (~0.125 
Gyr for A2199), leading to a cooling catastrophe in less than 
~3 Gyr. The authors found similar trends at different ini- 
tial magnetic field strengths and topologies. Starting with 
lower initial core densities and isothermal initial tempera- 
ture the cooling catastrophe was significantly delayed to ~7 
Gyr. The cooling catastrophe was also delayed when a cen- 
tral heating source was assumed via an idealized model of 
the heating luminosity within a radius of 20 kpc. The au- 
thors concluded (at least for A2199) that in the absence of 
AGN feedback and/or strong magnetic fields near the galaxy 
cluster core, thermal conduction alone cannot be responsi- 
ble for bal ancing the radiativ e losse s in the ICM, in agree- 
ment with IBogdanovic et all (|200gT ). Similar results about 



the reductio n and evolution of the ef fective heat flux were 
obtained by IRuszkowski fc Ohl (|2010l ) using the 3-D MHD 
code FLASH^~ 

There is a large body of numerical work on 
various aspects of AGN-ICM interact ions that has 
spanned purely hydrodynamic mode l s (e.g. Churazov et alj 
20011; iBriiggen fc Kaiserl l200ll. |20o3: IRevnolds et all |2002|; 
Omma et al. 2004 ; Zanni et alj 20051; IVernaleo fc Reynolds! 
20061). as well as MHD models (e.g. IRuszkowski et alj|2004l; 
Robinson etHl |2004| ; Ijones fc De Yound 120051 ; IDe Young 



20101 ) . Yet, a conclusive picture regarding the cooling flow 



problem based on n umerical simulations appears to remain 
elusive. For example. I Vernaleo fc Reynolds! ([20061 ) simulated 
direct injection of a supersonic jet into the ICM atmosphere. 
However, this model failed to maintain long-term thermal 
balance in the cluster, leading to hypotheses of heating by 
the dissipat ion of sound waves or the decay of g l obal ICM 
modes (e.g lOmma et al.ll2004D . Also. IDe YounJ (|201Ch ar- 
gue based on comparisons with experimental data that AGN 
outflows are strongly decelerated, generating turbulent sonic 
or subsonic flows due to their interaction with the surround- 
ing medium, which support such hypotheses. 

In this paper we present results from a series of 2-D nu- 
merical experiments aimed at assessing the driving physics 
and related sensitivities in a galaxy cluster. We have cho- 
sen to simulate A2199, a coolin g-core cluster that has been 
simulated by other authors (e.g.lZakamska fc Naravanll2003l ; 
IParrish et al1l2009l ; IRuszkowski fc Ohll2010l ). and for which 
observed data for the electron number density and temper- 
ature exist at a relatively close distance (R ~1 kpc) from 
the cluster center. The density at the core is an important 
parameter in MHD parameter studies that aim to assess 
mechanisms associated with the sustainment (or collapse) of 
the thermal balance in a cluster, since the radiative cooling 
rate depends strongly on it (oc n^). Because A2199 is a clus- 
ter that is often modeled by MHD simulations, we perform 
here a series of idealized numerical experiments in 2-D ax- 
isymmetric geometry to assess the sensitivity of the collapse 
rate on the near-core density, using two diff erent spatial pro- 
files that bounded the observational data ijjohnstone et al] 
[2002(1 . 

In 3-D MHD simulations it is common to prescribe ini- 
tial profiles that are in hydrostatic equilibrium. Often, the 
observational data guide solutions to the hydrostatic equa- 
tion and then these solutions are prescribed as the "ini- 
tial conditions" of the ICM. MHD simulations then pro- 
ceed to seek of mechanisms that maintain thermal balance 
in the cluster. In this study we perform a series of numer- 
ical experiments that do not begin with solutions in hy- 
drostatic equilibrium, to assess dynamical effects in the en- 
ergy balance of the system. Our simulations include a range 
of initial magnetic fiel d profiles (similar to th ose used by 
IParrish et al. 1 120091 and iBogd anovic et alj|2009l ) and quan- 
tify the self-consistent evolution of the ICM in 2-D planar 
geometry. 

The paper is organized as follows. Section [2] describes 
the computational methodology and governing laws used 
to perform the numerical simulations. Initial and boundary 
conditions for the cluster are also given in this section. Sec- 
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tion [3] presents our simulation results in 2-D axisymmetric 
(r,z) and in 2-D planar (x,y) physical domains. We summa- 
rize the results and provide our conclusions in Sections [4] 
and [5] respectively. 



2 METHOD 

The numerical simulations have been performed with the 
Multi block Arbitrar y Coordinate Hydromagnetic (MACH) 
code jPeterkirJll99gh . MACHz (where x — 2 or 3) are time- 
dependent, resistive, compressible MHD solvers of the Ar- 
bitrary Lagrangian Eulerian (ALE) variety that were de- 
veloped by the Center for Plasma Theory and Computa- 
tion at the Air Force Research Laboratory, and NumerEx 
Inc., to support applications with non-ideal MHD physics in 
complex geometries. MACH2 and MACH3 are the two-and- 
one-half-dimensional (2|-D) and 3-D versions of the solvers, 
respectively. Over the last three decades they have been up- 
graded significantly through the contributions of many sci- 
entists in the United States. Although the majority of the 
work with MACH today is in the area of pulsed power, their 
application base has expanded to multiple areas in engineer- 
ing. However, to the best of our knowledge, they have not 
yet been used to study problems in astrophysics. In order to 
introduce these codes to the astrophysical community, we 
provide a more detailed description of the codes and their 
applications in Appendix A. 

MACH2 is the version used in the present study. A 
feature of MACH2 that has been particularly useful in our 
studies is its modularity with the terms in the MHD equa- 
tions, which allowed us to perform instructive comparisons 
between the various simulation physics in the evolution of 
the ICM. The numerical experiments have been carried out 
using a reduced set of the full MHD conservation laws of 
MACH (see Appendix A). Simulations in only two dimen- 
sions allow a wide range of highly-resolved cases at reduced 
computational time, but at the expense of excluding pro- 
cesses in the evolution of the ICM that are inherently 3-D. 
Nevertheless, the range of physics cases examined has pro- 
vided us with a useful set of results, especially in regards to 
the "initial state" of the ICM and the potential importance 
of low-speed hydrodynamics near the center of the cluster. 
These 2-D results now also form the basis of our future sim- 
ulations with MACH3. 



2.1 Governing equations and numerical approach 

The ICM simulations have been performed using the follow- 
ing reduced set of the full MACH conservation laws, rep- 
resenting ideal MHD physics in the presence of anisotropic 
thermal conduction and radiative cooling: 



dp 
dt 



-V • (pv) 



o^j- = — pv ■ Vv — Vp + po 1 (V x B) x B + pg 



p— = -pv Ve — pVv — Vq- $ e H 
at 



8B 



— = V x (v x B). 



(1) 
(2) 
(3) 
(4) 



In equations Q-Q p is the mass density of the plasma fluid, 
v is the fluid velocity field, B is the magnetic induction field, 
e is the internal specific energy, q is the conductive heat 
flux, p is the thermal pressure and g is the gravitational 
acceleration. 

Analytic models or semi-empirical tabular equations 
of state (EOS) and transport coefficients may be used in 
MACH by accessing the sesame library (see Appenix A 
for more information about the library) during each compu- 
tational cycle. For the cluster simulations we have used an 
ideal-gas EOS with a specific heat ratio 7 = 5/3, 

Also, an analytic model of the anisotropic thermal con- 
ductivity coefficients has b een used, which incorporates the 
full Braginskii coefficients (Braginskii 1965) for the parallel 
(«W/) an d perpendicular (n e ±_) conductivities 

(6) 



K e //eseB + K eJ _ (I — e B e B 



where R e is the thermal conductivity tensor, &b is the unit 
vector in the direction of the magnetic field and I is the delta 
tensor; &B&B is the magnetic field unit dyad. The Braginskii 
coefficients for R e are provided in Appendix A. The conduc- 
tion heat flux is given by 



q = -K e ■ VT e 



(7) 



In the relevant range of temperatures (~1-10 keV) 
the plasma in the cluster cools primarily by thermal 
bremsstrahlung at a rate given by ( in units of J/m 3 /s) 
$ T b = 1.4 x W- 2a g b n e 2 T e 1/2 (e.g. iRvbicki fc Lightmanl 
1 19791) with (?;, denoting the velocity-averaged Gaunt fac- 
tor and n e is the electron number density. MACH2 uses 
a mean Planck opacity Xp (with dimensions of area cross- 
section/unit mass, so a photon mean- free-path is \/xp)- 
The opacity values are taken from the sesame library. In 
the optically-thin limit radiation energy is removed from a 
MACH2 computational element at a volumetric cooling rate 
of 



<t> e fl = aCiPXpTe 



(8) 



where a is Stefan's constant and ct is the speed of light in 
vacuum. A comparison of the cooling rates &Tb and & e R is 
provided in Appendix A. 

The gravitational acceleration is obtained from Pois- 
son's equation 



V • g » -ArxGpDM 



(9) 



where pdm is the mass density of the dark matter. New- 
ton's gravitational constant is denoted by G. The present 
study accounts only for the dark matter density, which is 
prescribed based on a Navarro-Frenk- White (NFW) profile 
with a soft core l|Navarro et al.|[l997h . The density profile 
may be expressed in non-dimensional form as follows: 



Pdm (r) 



(10) 



(f + f c ) {r + lf 

where p DM = 2nR a 3 p DM /M , f = R/R s and f c = R c /R s 
with R denoting the radius measured from the cluster origin 

2 http://tlweb.lanl.gov/doc/SESAME_3Ddatabase.1992.html 
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(r = 2 = or x = y — 0). The normalization mass is taken 
to be M = 3.8 x 10 14 M Q from lZakamska fc Naravanl (120031 ) 
and it has been kept fixed throughout all simulations pre- 
sented herein. The sta ndard scale radius for A2199, R a =390 
kpc, is taken also from lZakamska fe Naravanl (|2003l ) who in- 
ferred its value from observations of the outer plasma tem- 
perature. The value of the core radius R c is uncertain. For 
the ten clusters studied by Zakamska & Narayn the authors' 
"best-fit" solutions correspond to values of f c = and 1/20. 
For A 2199 they chose R c = 0, whereas iRuszkowski fc Oh] 
|2010l ) chose R c =20 kpc. We found that our main conclu- 
sions were not affected by this value and present our results 
for R c =20 kpc. The solution of equation ([9]) for the gravi- 
tational acceleration is spherically symmetric, g = g(R)e.R, 
where &r is the position unit vector defined relative to the 
cluster origin. In non-dimensional form the analytical solu- 
tion to equation ((9]) is given by: 

1 2r c - 1 , , 
9(r)= - -J-T^T, - + TT2 ln (^+ 1 ) 



f (f + 1) (f c - 1) 



■In 



(11) 



f 2 (fc - if 

where g = R s 2 g/2GMo- Based on estimated values of the 
NFW characteristic radii and normalized mass for A2199 we 
find Pbai ~ 1-3 x 10~ 24 g/cm 3 at R =0.95 kpc. At the same 
location the observed electron number density is 0.115 cm -3 
and is the highest value we have used in our initial condi- 
tions. This gives p ~ 0.2 x 10 -24 g/cm 3 for the plasma assum- 
ing that the molecular weight per electron is p e =1.18. The 
latter assumed that the hydrogen an d helium fractions are 
X = 0.7 and Y = 0.28, respectively (jZakamska fc Naravanl 
2003). Thus, our assumption in equation (J9]) about the dom- 
inance of the dark matter density in the determination of the 
gravitational acceleration appears acceptable. 

Equations (0-([4|) are solved in a fractional time-split 
manner. The Lagrangian time advance of the mass density, 
velocity and magnetic vector fields are computed implicitly 
on a computational domain that, unless otherwise noted, is 
discretized by a uniform mesh that consists of square com- 
putational elements (or "cells"). Simulations in both 2-D 
axisymmetric (r,z) and 2-D planar (x,y) coordinate systems 
have been performed but no mesh movement and/or adap- 
tation has been invoked. Spatial discretizations are obtained 
using finite volume differencing. The radiation cooling algo- 
rithm advances the internal energy explicitly, whereas ther- 
mal conduction is performed implicitly. Because we have as- 
sumed that the gravitational acceleration is dependent only 
on the dark matter density, g is prescribed using equation 
(|11[) and is held fixed throughout each simulation. 

2.2 Initial and boundary conditions 

2.2.1 Initial conditions for the plasma density and 
temperature 

Several combinations of initial conditions for the scalars and 
vectors have been used throughout this study, all associated 
with the A2199 galaxy cluster. Specifically, for the num- 
ber density and temperature of the plasma we employ two 
different profiles to illustrate the sensitivity of the ICM sim- 
ulations on the assumed initial conditions. The first profile 
is a fit of the Zakamska & Narayan "best" solution to an 



idealized set of four, time-independent, ordinary differential 
equations that account for thermal conduction and radia- 
tion cooling only, in hydrostatic equilibrium. The authors 
sought the solution th at best fits the Chandra observations 
l|johnstone et al.1l2002i ) by adjusting the thermal conduction 
heat flux factor /, defined here as 



/ = 



q ■ e R 



(12) 



where q sp = n e //VT e ■ ejj. The Spitzer-Braginskii thermal 
conductivity is K e /i = K sp = 3.16kB 2 n e T e T e /mc, where t c is 
the (classical) e-i collision time. The authors showed that for 
five out of the ten clusters studied a pure conduction model 
balancing radiative cooling, i.e. 



V • (/k sp VT e ) - $ eJ ? = 0, 



(13) 



was plausible if the factor / was reduced from unity by more 
than 50%. Because this best-fit solution yields an electron 
number density near the core which is about 35% lower than 
the observed value, we shall label it "LCD" (Low Core Den- 
sity). The LCD profiles are p lotted in Fig. [Q against Chandra 
data (|johnstone et al] |2002| . also in IZakamska fc Naravanl 
l2003h . 

Some MHD simulations in the literature have employed 
initial density profiles that underestimate the observed elec- 
tron nu mber density by fact ors of 5-6 near the cluster co re 
(e.g. see lParrish et al.l (|2009l ) and lRuszkowski fc Ohl l|201fj|) ). 
To assess the importance of the core conditions on the 
evolution of the ICM we constructed a second set of ini- 
tial conditions in a ma nner similar to that followed by 
IRuszkowski fc Ohl ([2010) , who assumed hydrostatic equilib- 
rium and a pol ytropic EOS with a pr escribed profile for the 
entropy S (e.g. Cavagnolo et al.ll2009l ). We write the entropy 
as a function of radius in non-dimensional form: 



S(f) 



07-1 



1 + fif/c, 



(14) 



where the barred quantities are defined as S = S/S c , 
f = Tc/Tc, p = p/pc, 13 = (T 20 o/Tc)(p c /P20o) 7 " 1 and 
c = -R200/ 'Rs is the concentration parameter. Here, variables 
carrying subscript "c" represent values at the cluster center 
and those with a subscript "200" denote values at a radius 
where the mean density of the dark matter is 200 times the 
critical density of the universe at the redshift of the cluster 

pcrit(z). 

The hydrostatic equation is given by 



Vp = pg. 



(15) 



Assuming spherical symmetry, equation (|15|l may be writ- 
ten in terms of the non-dimensional density and entropy as 
follows: 



dp 
df 



75 



Cg 



df J 



(16) 



and, given S(f), yields the density as a function of ra- 
dius. The gravitational acceleration is taken from equa- 
tion pip . The constant on the right of equation (|15p is 
C = 2GMo/j,m u / RsksTc (/i denotes the mean molecu- 
lar weight, m u is the atomic mass unit and ks is Boltz- 
mann's constant). Equation (|15[1 has been solved numeri- 
cally to produce a second set of cluster profiles by assum- 
ing the following parameters: T c =1.6 keV, T200 =4 keV, 



2-D MHD Simulations of Galaxy Clusters 5 



10- 



> 



1 - 



T e Chandra data 
n e Chandra data 
T e HCD model 
T e LCD model 
n e HCD model 
n e LCD model 




f 10 



10 



10 



r 10" 



0.1 



10 

R (kpc) 



100 



Figur e 1. Chandra da ta l|johnstone et alj l200l also in 
Zaka mska &: Nar avan 2003) and initial profiles used in the A2199 
galaxy cluster simulations. 



p c = 23.5 x 10 - 26 g/cm' 3 , p 2 oo = 0.054 x 10~ 26 g / cm~ 3 , 
c =4.378, 7 =5/3 and /1 =0.62. The profiles are plotted 
in Fig. Q] and shall be denoted "HCD" (High Core Den- 
sity) hereinafter. Although neither the LCD nor the HCD 
profiles represent accurately the density observations, they 
were implemented in this manner deliberately as limiting so- 
lutions that bounded the full range of the data, within the 
constraints of the hydrostatic equilibrium condition, while 
illustrating the sensitivity of the simulation results on the 
density. 



2.2.2 Initial conditions for the magnetic field 

An extensive body of work exists on the magnetic field 
structure in galaxy clusters covering both cooling-core and 
non-cooling-core clusters. Faraday Rotation Measures to- 
wards extended radio sources in and behind galaxy clusters 
have demonstrated the turbulent natu r e of the magnetic 
field there (e.g. see iGovoni et al.l |2010|. IVacca et ail 20101. 
Guidetti etaiTl2010l. IrBonafede et al.1 l20ll iGuidetti et al ' 



20081 lLaing et al.ll2006l . IVogt fc EnMrJ 12003 , iMurgia et al 
20041 ). In general it has been found that a power- 
law distribution of their fluctuation spectrum and some 
radially-declining profile of their average strength rep- 
resent well these magnetic field structures. Nearly all 
hydrodynamical simulations within a cosmological con- 
text pr edict around 10% of turbulence within the ICM 
e.g. seelBrvan fc Normanlll998l.llnogamov fc Sunvaevl2005 



or former groups of galaxies that move through the cluster 
with speeds that are typcially on the order of a thousand 
km/s, permanently perturbing the underlying potential. It 
has been shown that such gravitational perturbations alone 
can pro duce significant motions within the core of clusters 
(e.g. seelXscasibar fc Markevitchll2006l . iRoediger et alj20ld . 
IZuHone et alj|2010h . Some very small fraction of these sub- 
structures contains also gaseous atmospheres or disks capa- 
ble of steering the ICM directly. 

To assess the impact of initial magnetic field topolo- 
gies on the evolution of the effective thermal conduction 
in the ICM we employed here the following three initial 
magnetic field configurations: (1) a turbulent field (using 
random number deviates), (2) a purely rotational field and 
(3) a purely radial field. The intent of the last two highly- 
idealized configurations was to assess the sensitivity of the 
anisotropic heat flux reductions, and have served as the two 
limiting cases of our parameter studies on magnetic fields. 
In all three cases the out-of-plane component of the mag- 
netic field, Be in the 2-D axisymmetric cases and B z in the 
planar cases, was set to zero. 

The turbulent fie ld was constructed bas ed on the 
approach followed by IRuszkowski et al.l ((2007j) (see also 
iRuszkowski fc Ohll2010l ). First, we set up a random field 
in fc-space, using independent normal random deviates 
^(0,1) (with and 1 being the values of the mean and 
standard deviation of the random distribution, respectively) 
for the real and imaginary components: 



B x (k) = [Re{B x ),Im{B x )\ 
B y (k) = [Re(B y ),Im(B y )] 



[NiB,N a B] 
■ [N 3 B,N 4 B] 



(17) 

where the wavenumber magnitude is fc = (k x 2 + ky 2 ) 1 ^ 2 
and subscript "a" is a number used to identify different ran- 
dom deviates. Ruszkows ki et all l|2007l ) proposed the follow- 
ing dependence of the magnetic field amplitude on the wave 
number \k\: 



B(k) 



-11/6 



exp [— (k/ko) 4 ] exp (— ki/k) . 



(18) 



where ko and fci control the exponential cutoff terms in the 
magnetic energy spectra. The former controls the spectra at 
the large wave numbers and is given by ko = 2n/Xo with 
Ao oc 43ft -1 and h denoting the normalized Hubble param- 
eter (0.5 — 0.8). The second, ki, controls the cutoff of the 
spectra at the small wavelengths; the case we present in this 
paper uses ko = 0.095 and fci = 0.05. Once the complex 
vector field in fc-space has been generated, we perform vec- 
tor divergence cleaning. Using conventional Einstein tensor 
summation notation where Su is the Kronecker delta, this 



Gardini et al.ll2004 IVazza et alj|2006l . Ilapichino et aliboosl. operation may be expressed by (e.g. see| Balsara||l998| ): 



Paul et alj |2010h. and are in agreemnt with observations 



(e.g. see lSchuecker et al.ll2004 IChurazov et al.ll2004l ). Con- B % {k) = I 8 



sequently, since the magnetic field is frozen within such tur- 
bulent ICM, cosmological MHD simulations predict that the 
magnetic field h as a turbulent spe c trum (specific examples 
can be found in iDolag et al. 2002, iBriiggen fc Hoeftl I2005T ) 
with no exception when studying a sample of clusters or 
during the time evolution of a single cluster. The reason 
for this is that clusters forming in a cosmological context 
contain accreting matter, part of which forms by small- 
scale-structure mergers. Therefore, clusters may be popu- 
lated by thousands of sub-structures in the form of galaxies 



(19) 



Finally, we perform (complex) inverse Fourier trans- 
formation in 2-D to obtain the r-space components 
[B x (r),B y (r)]. 

Because the precise distribution of magnetic fields in 
galaxy clusters is uncertain, there is no strong evidence in 
favor of the assumed spectrum in equation (| 18f) . Therefore, 
we extend the range of our numerical experiments to include 
the following two (extreme) cases: a rotational profile and 
a radial profile as given by equations (|20l) and (|21l) . respec- 



6 I. G. Mikellides, K. Tassis and H. W. Yorke 



500- 
400- 
300- 



I ■ :-■ ■ 



-500 -^ffl! 




\ <U|>c> 



; (kpc) 



0.01- 

(V 

4- 



0.00 M 




■ Turbulent 

■ Radial & Rotational 



500- 
400- 



i 1 1 1 1 1 1 1 1 1 r 

50 100 150 200 250 300 350 400 450 500 • 
x (kpc) 







-3-»>V' '"■'-'"' ' 




- -17'' *W _ 

r r - ^' ' ' - ■ . 












p i ti-|VCfi i n i i'i 1 


iiiLrrmv^ i i ,ir, --M^v ( |iiri.^Vi i' 



son -joo -300 



; f,i',v 
[,9x10 

6 5x10' 
22x10 

7 5x10" 
2,5x10" 
£E.6xlO 
2,9x10" 
I.OxlO 



Figure 2. Initial magnetic field profiles used in the 2-D planar simulations. The x-y plots show contours of magnetic pressure (in Pa) 
overlaid by magnetic field unit vectors. The magnitudes of the magnetic field for the three cases used in the simulations - rotational, 
radial and turbulent - are compared in the lower left plot along x at y=0. In all planar (x,y) simulations the computational mesh was 
uniform consisting of 128 2 square cells. 



tively, where Ro is a scaling radius. 
B *(*,V) = -Bj^- 2 ,B y (x,y) = +B-^±- 2 



ar + y 2 - 
xRo 



ar + y 1 



(20) 



(21) 



B x {x,y) B x * +y i, aP+y* 

These two additional cases have been used exclusively in 
the 2-D planar simulations and the results are presented in 
Section [3.21 For both cases the out-of-plane component B z 
was set to zero. The initial magnetic field profiles for all 
three cases are illustrated in Fig. [5] for a peak magnetic field 
value of 1 fiG. Hereinafter we shall denote the maximum 
strength of the magnetic field in our simulations by B. 



2.2.3 Boundary conditions 

Boundary conditions in MACH2 are defined by employing 
"ghost" cells and vertices. The simulations in planar ge- 
ometry were performed using a computational domain that 
spanned — L < x < L and — L < y < L with L —500 kpc. 



At all outer boundaries in this geometry we have imposed 
a Dirichlet condition for the temperature and the specified 
value is set equal to the initial value as given in Fig. [T] Ra- 
diative cooling is allowed at these boundaries. For each time 
step the velocity at the outer boundaries is set to zero ( "no 
slip"), the density is set to the initial value (QJ, an d the 
magnetic field is copied from the real cell adjacent to the 
boundary to the ghost cell. 



For the axisymmetric (r,z) simulations the computa- 
tional region spans < r < L and < z < L. At the 
boundaries (L, z) and (r, L) the conditions are the same as 
those at the outer boundaries of the planar geometry. At 
(0, z) and (r, 0) we impose thermally insulating conditions. 
Only one case in this computational domain employs a mag- 
netic field and it is kept frozen everywhere throughout the 
simulation. 
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Table 1. Numerical experiments in 2-D axisymmetric (r,z) ge- 
ometry. Only thermal conduction and radiation were included in 
the simulation physics (hydrodynamics were not allowed in these 
cases). The heat flux factor / was either specified or allowed to 
evolve freely. 





Initial Conditions 


Heat Flux 
Factor, / 


n e & T e 


B(pG) [Profile] 


LCD 


0.0 [N/A] 


0.4 


LCD 


0.0 [N/A] 


0.5 


LCD 


0.0 [N/A] 


0.6 


LCD 


0.0 [N/A] 


1.0 


LCD 


1.0 X 10~ 12 [turbulent & frozen] 


free 


LCD 


0.1 [turbulent & frozen] 


free 


LCD 


1.0 [turbulent & frozen] 


free 


HCD 


0.0 [N/A] 


0.5 


HCD 


0.0 [N/A] 


0.6 



Table 2. Numerical experiments in 2-D planar (x,y) geometry. 
The heat flux factor / was allowed to evolve freely. (Simulation 
physics nomenclature: "C"= Thermal conduction, "R"= Radia- 
tion cooling, "H"= Hydrodynamics, "M" = Magnetohydrodynam- 
ics.) 





Initial Conditions 


Simulation 
Physics 


n e & T e 


BQjG) [Profile] 


LCD 


0.0 [N/A] 


CR 


LCD 


0.0 [N/A] 


CRH 


LCD 


1.0 [turbulent & frozen] 


CR 


LCD 


1.0 [rotational] 


CRM 


LCD 


1.0 [radial] 


CRM 


LCD 


1.0 [turbulent] 


CRM 


LCD 


1.0 X 10~ 3 [turbulent] 


CRM 


HCD 


1.0 [turbulent] 


CRM 


HCD 


1.0 X 10 -3 [turbulent] 


CRM 



3 NUMERICAL EXPERIMENTS 

A series of numerical experiments were performed in 2-D 
axisymmetric (r,z) and 2-D planar (x,y) geometry with the 
two initial plasma profiles described in Section 12.2.11 LCD 
and HCD, and three initial magnetic field profiles described 
in Section [2. 2. 2 1 for a wide range of combinations and simu- 
lation physics. Hereinafter, it will be implied that when we 
present simulation results without a magnetic field, ther- 
mal conduction has been computed isotropically, whereas 
for simulations that include a magnetic field thermal con- 
duction has taken into account the anisotropic effects (see 
equations [S] and [7} . 

Table [T] lists the cases carried out in 2-D axisymmet- 
ric geometry. A goal of the simulations in this geometry 
was to assess the sensitivity of the global heat flux factor / 
on the plasma density near the A2199 core under idealized 
simulation physics. Specifically, no hydrodynamics were per- 
mitted in any of these simulations and the magnetic field, 
when used, was held frozen. Table [5] outlines our cases in 
2-D planar geometry. In all of these cases no control of the 
thermal conductivity was exercised and the anisotropic heat 
flux was computed self-consistently. The focus of these sim- 
ulations was the early dynamics in the near-core region of 
the cluster and the results are presented in Section \3. 21 



3.1 2-D simulations in axisymmetric (r,z) 
geometry 

The main intent of the simulations in 2-D axisymmetric 
geometry was to conduct a preliminary assessment of the 
sensitivity of the ICM solution to the prescribed initial 
conditions before invoking full MHD simulations. Due to 
the dependence of the cooling rate on the square of the 
electron density, $ e _R oc n 2 \/Tl, we expect the solution 
will depend strongly on the near-core gas density (e.g. see 
IConrov fc Ostrikerl jjooj ) . In fact, if one assumes that the 
dominant competing mechanisms establishing the thermal 
balance in the ICM are thermal conduction and radiation 
only, and all other mechanisms are captured by a single ef- 
fective reduction of the thermal conductivity (represented by 
the quantity /), we can estimate of this sensitivity from the 
steady-state energy equation [T3] by computing directly the 
heat flux factor / of the cluster. The value of the ICM den- 
sity in the cluster core (R ~0.95 kpc), based on the Chandra 
observations, is n e « 0.115 cm -3 , whereas the LCD model 
at the same location gives n e w 0.075cm -3 . Since the tem- 
perature predicted by the LCD and HCD models is approx- 
imately the same at this location, we use the LCD model to 
estimate V-(K sp VT e ). By approximating the core region as a 
small finite cylindrical volume with dimensions Ar = Az « 
L/128=3.906 kpc we find: 



(i) fas 0.8 when n e 

(ii) fas 0.5 when n e 



0.115 cm~ 
i 0.075 cm" 



which illustrates to first order the sensitivity of idealized 
solutions on the near-core plasma density. 

The significance of this sensitivity may be better appre- 
ciated by comparing the following relevant time scales. First, 



tion of a region with a temperature gradient 
may be estimated by setting 



3 , 9T« 
2 nekB —t 



AT e /AR, 



(22) 



yielding 
3 



, AR 2 n e t 2 
n e k B — oc 



(23) 



(with £ ~ AR). Next, the characteristic time rj> in which the 
plasma radiates away all of its internal energy is obtained in 
a similar manner yielding 



3 n e k B T e 
2 actpxpTe 



-a/2 



(24) 



If the mechanism(s) in the ICM produce an effective 
reduction of the thermal conductivity fm 0.5, then the con- 
duction time would be more than double the radiation time, 
assuming the plasma density near the core is as high as 
the observed value. In other words, without any additional 
heating source(s), these estimates suggest that thermal con- 
duction from the outer regions of the cluster cannot avert 
runaway cooling near the core. 

Indeed, 2-D axisymmetric simulations with MACH2 
that incorporated only thermal conduction and radiation 
(MHD was turned off in this series of calculations) show that 
thermal balance is achieved at a relatively low value of / (be- 
tween 0.5-0.6) for the LCD profile, whereas the HCD profile 
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Figure 3. Results from 2-D axisymmetric simulations with ther- 
mal conduction and radiation only. The solutions are plotted for 
different values of the heat flux factor /, at R =20 kpc. Top: Re- 
sults using the LCD initial conditions for 64 2 and 128 2 uniform 
(U) computational meshes. The 64 2 non-uniform (NU) mesh used 
quadratic spacing with highest resolution of 0.25 kpc at the clus- 
ter core. Bottom: Comparisons of solutions for the LCD and HCD 
initial conditions at f=0.5 and 0.6 (it is noted that differences in 
the two HCD solutions are negligible). 




Figure 4. Results from 2-D axisymmetric simulations with 
anisotropic thermal conduction, radiation, and a turbulent mag- 
netic field. The magnetic field was frozen throughout the simula- 
tion. Top: Temperature contours overlaid by magnetic field unit 
vectors at t=0. Bottom: Evolution of the average temperature at 
R =20 kpc for different values of the peak magnetic field ampli- 
tude allowed. It is noted that differences between the solutions 
for the 0.1 and 1 fiG are negligible. 



leads to the collapse of the core for the same values of /. Fig- 
ure[3](top) shows the evolution of the computed temperature 
at R =20 kpc for different values of /, when LCD initial pro- 
files were used. Mesh sensitivity calculations confirmed that 
a 64 2 uniform mesh provided sufficient resolution. Figure [3] 
(bottom) compares the LCD solutions to those obtained us- 
ing the HCD profiles for / = 0.5 and 0.6, showing collapse 
of the core in less than 1 Gyr. Thus, determining a global 
value for / that would sustain thermal balance in the cluster 
may be in general instructive; quantitatively, however, it is 
of less significance due to the sensitivity of / on the central 
density. 

Due to the extremely high Hall parameters in the ICM 
(fl e > 10 9 in A2199) the potential impact of turbulent mag- 
netic fields in the transport of heat from the outer (hot) re- 
gions of the cluster has been widely recognized. Although the 
evolution an d impact of MHD instabilities, specifically o f the 
HBI (e.g. see lQuataertll2008l : |Parrish fe QuataertfeOOSt) and 
of the magneto-thermal instability (MTI) (e.g. see Balbusl 
l2000l : |Parrish" et al.l [2 008l ) , remain an ongoing topic of study, 



we performed here an idealized numerical experiment with a 
turbulent initial magnetic field (see Section [2. 2. 2[) that was 
kept frozen throughout the simulation. The intent was to 
emulate a scenario in which continuous galaxy motions sus- 
tain a nearly-isotropic turbulent field. T he simulation case 
was m otivated in part by the work of iRuszkowski fc Oh] 
(2010), who suggest it is possible that turbulent velocity 
fields in the cluster may drive (or stir) turbulence in the 
magnetic field in a manner that allows for sufficiently rapid 
transport of heat from the outer regions of the ICM to sus- 
tain thermal balance. 

Figure|4](top) shows contours of the initial temperature 
overlaid by unit vectors of the initial magnetic field as im- 
plemented in our 2-D axisymmetric simulation. For this case 
the peak magnitude of the magnetic field did not exceed 1 
juG. The evolution of the computed temperature at R — R c 
is shown in Fig. [3] (bottom). Also shown for comparison is 
the case with no magnetic field. For all cases in this series 
of simulations the thermal conductivity was computed self- 
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consistently, accounting fully for the anisotropy in the heat 
flux, but magnetohydrodynamics were excluded. We found 
that thermal conduction alone is unable to avert a cooling 
catastrophe, even when the magnetic field is completely tan- 
gled. The collapse of the core for the 1-/^G case was found to 
be insensitive to the value of peak magnetic field, which is 
expected due to the high values of f2 e . In this arrangement a 
change in the collapse rate was obtained only when the peak 
magnetic field magnitude was lowered by several orders of 

magnitude (where Q e —¥ 1). 

By comparison, iRuszkowski fc Obi l|2010h found in their 
3-D simulations of A2199 that a cooling catastrophe could 
be averted. The authors applied a driving force that resulted 
not only in keeping the magnetic field tangled but also in 
turbulent flows that introduced additional physical effects 
such as mixing, which we do not account for in our 2-D ax- 
isymmetric simulation. We also recognize that the real topol- 
ogy of the magnetic field is unknown and may differ signif- 
icantly from the idealized model described in Section \2. 2. 2 1 
Nevertheless, in view of the sensitivity of thermal balance 
calculations on the plasma density it appears that hypothe- 
ses about "stirring" or other dynamical activity in the ICM, 
proposed as mechanisms that could sustain thermal balance 
in the cluster, would only be strengthened if such mecha- 
nisms sustained balanc e in the cluster for a wid e range of 
possible core densities. Ruszkow ski fc "Ohl l|2010h employed 
an initial temperature profile that was in close agreement 
with the Chandra data and a plasma density profile with 
values almost six times lower than the observed data near 
the cluster core. 



3.2 2-D simulations in planar (x,y) geometry and 
in hydrostatic non-equilibrium 

In this section we perform an idealized study in two dimen- 
sions near the core of A2199, focusing on po ssible mecha- 
nism s that may overcome HBI -driven collapse (jParrish et al.l 
120091 ; IRuszkowski fc Ohll2010l '). In the 3-D MHD simulation 
cases of iParrish et all (|2009h and IRuszkowski fc Ohl (|2010l ) 
some initial ICM states did not achieve thermal balance but 
all cases imposed hydrostatic equilibrium at t=0. Here, we 
deliberately violate this condition only as a means to gener- 
ate and study near-core flow dynamics (<20 kpc) that may 
possibly be associated with AGN activity. 

We deviate from hydrostatic equilibrium by prescribing 
the spherically-symmetric initial conditions LCD and HCD 
(Fig. [Q in a 2-D planar geometry. We note that although 
both these profiles correspond to hydrostatic equilibrium in 
three dimensions, they do not preserve such equilibria in two 
dimensions, thereby producing subsonic wave motion in the 
plasma. Here we focus only on the early dynamics of the 
MHD flow field, t < 2.5 Gyr, near the cluster center (R <20 
kpc), since several simulations of A2199 have predicted a 
cooling catastrophe due to the HBI within this time for a 
wide range of magnetic field strengths and profiles. Assum- 
ing a weak magnetic field initially in the direction of the 

ravitational for ce, the HBI has a characteristic growth rate 

Quataertll2008h 
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Figure 5. Comparisons of the average temperature at R =20 
kpc, computed in planar geometry for cases with the same (LCD) 
profiles but different simulation physics (see Tablc[2]for simulation 
physics nomenclature). For all cases with a magnetic field the 
maximum strength was 1 /iG. 



~ (0.125 Gyr) -1 . Thus the instability is allowed approxi- 
mately 20 growth times in ~2.5 Gyr. 

We find in the 2-D planar simulations with a turbulent 
initial magnetic field (see Fig. [2] bottom right) that the av- 
erage temperature at R =20 kpc is sustained and a cooling 
catastrophe does not occur, as shown in Fig. [5] For com- 
parison the solutions with the same simulation physics are 
also shown in Fig. [S] for the following cases: (1) radial and 
(2) rotational initial magnetic fields, (3) turbulent initial 
magnetic field profile but held frozen throughout the sim- 
ulation (a case with no MHD but one that accounted for 
anisotropic thermal conduction and radiative cooling), (4) 
isotropic thermal conduction and radiation without MHD, 
and (4) isotropic thermal conduction, radiation and hydro- 
dynamics but no magnetic field. Similar to the cases of 
isotropic conduction no catastrophic cooling is found in the 
simulations with large-scale flows (induced by starting out 
of hydrostatic equilibrium), independent of the initial con- 
figuration of the magnetic field and despite a significant re- 
duction of the heat flux factor. Figure[U](left) plots the evo- 
lution of the cosine of the angle between the magnetic field 
and radial position vectors, (\e.B ■ e-n\) and the heat flux fac- 
tor (f) in the case of the turbulent initial magnetic field and 
the LCD profile. Unless otherwise noted all quantities have 
been averaged over the area of a ring of radius 20 kpc and 
maximum thickness of A, which denotes the side length of a 
single computational cell (recall our mesh consists of square 
cells). Also plotted in Fig. [6] (left) are the absolute radial 
velocity ([1; • en\) and Mach number (\v ■ e.n\) / (c s ) where c s 
is the adiabatic acoustic speed. 

The evolution of the velocity field as shown in Fig. [6] 
(left) defines a characteristic interval of time, 0.2 < t < 1.3 
Gyr. The interval is bounded by two characteristic times. 
The shortest time is associated with the hydrostatic imbal- 
ance that is imposed initially in the system, producing a sub- 
sonic expansion with a characteristic speed of ~100 km/sec 
((M) ~0.1). The time it takes for the expansion wave to 
reach a radius of 20 kpc is ~20 kpc/100 km/sec=0.2 Gyr. 
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Figure 6. Results from 2-D planar simulations with CRM simulation physics starting with the LCD initial condition for the density 
and temperature and a turbulent initial magnetic field with maximum magnitude of 1 /iG. Left: Evolution of the cosine of the angle 
between the magnetic field and radial position vectors, {\eg ■ &r\), heat flux factor (f), absolute radial velocity (|«-efl|) and Mach 
number (\v ■ en\) / (c s ) , averaged over the area of a disc of radius R=20 kpc (±A), Right: Temperature contours (in eV) overlaid by 
magnetic field unit vectors at t=1.27 Gyr showing the orientation of the field to be largely rotational near the cluster core. 




x (fcpe) x (hue) x (kne) 

Figure 7. Evolution of the temperature and magnetic (unit vector) field for the case of the LCD initial profile and a 1 fiG turbulent 
initial magnetic field. The calculations employed CRM simulation physics in planar geometry. 



As shown in Fig. [B] (left), upon passage of this first wave the 
rarefaction flow at R =20 kpc moves at an approximately 
constant speed, ~10 km/s ((M) ~0.01) for about another 
1.1 Gyr. Such "winds" are, in fact, observed in luminous 
galax ies (and are comm only labeled broad-line regions, e.g. 
sec lOstriker et "all l2010l l . In the ensuing discussions we fo- 
cus in this "quiet" rarefaction period 0.2 < t < 1.3 Gyr as 
an emulator of low-speed flow activity (possibly associated 
with an AGN). Figure[7]plots 2-D temperature contours and 
unit vectors of the magnetic field at t=0, 0.16, and 0.32 Gyr 
showing the evolution of the expansion. 

The two-way transit time for sound waves generated 
near the core is ~ 2L/c s = 1.3 Gyr. At about this time the 
heat flux factor and cosine of the magnetic field angle at 
R =20 kpc have reached their minimum values of 0.17 and 
0.35 (angle=70 deg), respectively (note that the angle in 
fact decreases within t ~0.2 Gyr to about ~45 deg from 
the initial value of ~60 deg due to the expanding flow). 
Figure [6] (right) shows the re-orientation of the magnetic 
field lines due to the HBI within ~100 kpc at 1.27 Gyr, 
which i s in qualitative agreem ent with the 3-D simulation re - 
sults of iParrish et~aT1 l|2009l ) and iRuszkowski fc Ohl ( |2010h . 



Also noted is the qualitative similarity of our solution for 
the heat flux factor in F ig. [6] (left) with that obt ained by 
IParrish et al. I (|2009T i and IRuszkowski fc Ohl |2010l ) in three 
dimensions before saturation (to the value of 0.07). At about 
1.3 Gyr the velocity at R =20 kpc begins to increase as the 
first expansion wave is reflected off the outer boundaries of 
our physical domain and reaches the near-core region. Be- 
yond this time the magnetic field orientation with respect 
to the radial position vector, the heat flux factor, and the 
temperature are driven by the hydrodynamics of the flow. 

A closer look at the various terms in the energy equa- 
tion shows the diminishing importance of thermal conduc- 
tion near the core in the early energy balance of the system, 
with the convective and radiative contributions almost fully 
driving the system by ~1.4 Gyrs. Figure [8] compares the 
power densities for radiative cooling $ e j?,, thermal conduc- 
tion V • q and the convective terms -(pv ■ Ve + pV • v), at 
three different times along the y =0 axis. In all cases we 
found that pV ■ v is the dominant term over pv ■ Ve. Also, 
because the initial conditions imposed zero velocity, the plot 
at t — 3.2 x 10 -5 Gyr simply reflects that the evolution of the 
internal energy of the system is initially driven by thermal 
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Figure 8. Comparison of the contributions to the ICM energy balance from the convection and pdV work, radiation, and anisotropic 
thermal conduction. The comparisons are for the case of the LCD initial profile in planar geometry, employing CRM simulation physics 
and a 1 fiG turbulent magnetic field. The slice plots are along x at y=0. The computational mesh consisted of 128 2 square elements 
yielding a minimum size of 7.8 kpc for each element. 
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Figure 9. Evolution of (\e.B ■ (/), (\v ■ &r\) and (\v ■ ea\) /{c B ), averaged over the area of a disc of radius R=20 kpc (±A). The 

cases employed CRM simulation physics in planar geometry. Left: Results starting with the LCD initial condition for the density and 
temperature and a rotational initial magnetic field with maximum magnitude of 1 /iG. The velocity flow field here fluctuates between 
radially inward and outward directions after the first expansion at t~0.1 Gyr (note that we plot the absolute value of the radial velocity). 
Right: Results starting with the HCD initial conditions for the density and temperature and a turbulent initial magnetic field with 
maximum magnitude of 1 fiG. The flow field is directed inward after the first expansion at t~0.1 Gyr (see also Fig. Illll . 
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Figure 10. Evolution of the temperature (in eV) and magnetic (unit vector) field for the case of the LCD initial profile and a rotational 
initial magnetic field with peak magnitude of a f fj,G. CRM simulation physics were used in planar geometry. 
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conduction and radiation, but by ~1.4 Gyrs temperature 
gradients near the core are found to be almost completely 
smoothed out. 

The diminishing effect of anisotropic thermal conduc- 
tion is even more evident in the case of a rotational initial 
magnetic field, which is found to reduce on- average the heat 
flux to the core within the time interval of interest by more 
than one order of magnitude compared to the case of a tur- 
bulent initial field (see Fig. [6] (left)). The evolution of the 
pertinent average parameters is shown Fig. [9] (left). The ro- 
tational magnetic field remains approximately orthogonal to 
the radial unit vector; the angle between them does not fall 
below ~78 deg. Yet, this has little effect on the evolution 
of the temperature during this time. This may seem coun- 
terintuitive at first since the characteristic times associated 
with the smoothing of temperature gradients in the 1CM 
by thermal conduction are larger in the rotational-field case 
as the effective heat flux factor is significantly lower com- 
pared to the turbulent field (equation 1251 gives r q ~f2 Gyr 
for /=0.02). However, thermal conduction does not drive 
the temperature in these cases. The characteristic time for 
temperature changes induced by the pressure work on the 
fluid may be estimated by 



(26) 



(27) 



3 U Y7 

-n e k B — — = -pV ■ v 
2 at 

yielding 

~ 3 ATe 1 
Tp ~ 2~tT~Av' 

For the induced conditions in our numerical experiments 
with the LCD profiles r p rsO.13 Gyr (<< r q ). Figure [TO] de- 
picts the evolution of the temperature showing the smooth- 
ing of the temperature gradients despite the preservation 
of a highly-rotational magnetic field within the time inter- 
val of interest. Thus, under the induced conditions pV ■ v 
dominates over thermal conduction. 

The 2-D axisymmetric simulations in Section [3.f I with 
only radiation and thermal conduction demonstrate a strong 
sensitivity of the solution on the assumed initial core den- 
sity, with the HCD profile yielding a collapse of the core at 
values of /for which the LCD profile yielded no cooling. The 
results of the cosine of the magnetic field angle, heat flux fac- 
tor, radial velocity and Mach number are depicted in Fig. 01 
(right) for the HCD profiles. The hydrodynamics of the flow 
here are significantly different compared to the LCD case. 
The induced hydrostatic non-equilibrium state for the ICM 
encompasses higher energy density in the system but at a re- 
duced pressure gradient. The resulting expansion is weaker 
compared to the LCD case (fO km/s wave front speed in 
HCD versus fOO km/s in LCD) and the gravitational force 
is able to turn the flow at R —20 kpc inward relatively fast, 
within about O.f-0.2 Gyr, and the plasma beyond this time 
continues to flow radially inward. This is in contrast to the 
LCD case with a rotational magnetic field where the veloc- 
ity field fluctuates between inward and outward directions 
(note that we plot the absolute value of the radial velocity). 
The temperature contours at t=0.f6 Gyr are shown in Fig. 
If f I Thus, the response of the temperature in Fig. If 21 and of 
the non-dimensional quantities in Fig. [9] (right) responded 
considerably faster than in the LCD case. The solutions with 
a peak strength for the initial magnetic field of f nG show 
only minor differences with the 1-fiG cases. 
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Figure 11. Temperature (in eV) and velocity (unit vector) field 
for the case of the HCD initial profile and a turbulent initial mag- 
netic field with peak magnitude of a f fiG and CRM simulation 
physics in planar geometry. 
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Figure 12. Comparisons between the average temperature for 
the LCD and HCD profiles starting with a turbulent initial mag- 
netic field with maximum magnitude of 1 nG and 1 fiG. All cases 
employed CRM simulation physics in planar geometry. 



4 SUMMARY AND DISCUSSION 

The works of IZakamska fc Naravanl l|2003h . 



IConrov fe Ostrikerl (|2008l ) and iParrish et all (|2009l ) have 
been particularly influential in formulating the different 
cases presented in this paper. Our 2-D axisymmetric 
si mulations reproduce cl o sely t he effective heat flux factor 
of IZakamska fc Naravanl l|2003h under similar assumptions 
for A2199, when thermal conduction and (bremsstrahlung) 
radiation are the only physics allowed in the energy balance 
of the cluster, and when the electron number density and 
temperature profiles closely resemble the authors' "best-fit" 
solution; we have termed this solution as the "LCD profile." 

Because the LCD profile underestimates the Chandra 
density observations near the core by approximately 35% 
(albeit in excellent agreement with the data beyond ~10 
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kpc), we performed a second set of simulations with an ini- 
tial profile that was in agreement with the observed electron 
density near the cluster center but overestimated the data 
at larger distances. Differences in the temperature profiles 
between the two cases were less than 5% for R<200 kpc. Be- 
cause the ratio of characteristic time for thermal conduction 
over the time for radiative cooling is strongly dependent on 
the electron density, t 9 /t$ oc n e 2 (assuming everyting else 
fixed), the goal of these simulations was to determine what 
effect such an increase would have on the thermal balance in 
the cluster. By comparion to the LCD case we found that the 
ICM collapsed catastrophically in the HCD case for the same 
values of the heat flux factor /. Our conclusions from this se- 
ries of idealized simulations regarding the importance of the 
core density in the ev olution of the ICM are in agreement 
with the conclusions of lConrov fc Ostrikerl ((2008) , who illus- 
trated the dependence of the ICM collapse rate on the value 
of the plasma density at the center. A turbulent initial mag- 
netic field that was held fixed throughout the simulation also 
led to the collapse of the ICM for both initial density profiles 
(LCD and HCD). In this series of simulations the anisotropy 
of the conduction he at flux was accou nted for by invoking 
the approximations of lBraginskiil 0965) for the thermal con- 
ductivity coefficients. T he magnetic field was ge nerated by 
a method proposed by iRuszkowski fc Ohl {2010 ') , and was 
held frozen to emulate an extreme case in which galaxy mo- 
tion maintains the turbulent field in the absence of any flow 
dynamics in the ICM. 

Since it remains largely infeasible to perform fully- 
consistent cosmological simulations of clusters that account 
for all pertinent physics and scale lengths, most often in 
MHD simulations the thermodynamic state and magnetic 
field topology are prescribed as the "initial conditions" of 
the ICM cluster and the simulations then proceed to seek 
of mechanisms that maintain thermal balance in the clus- 
ter. Such simulations do not always begin with an ICM in 
thermal balance but all studies we have reviewed imposed 
hydrostatic equilibrium at t=0. Referring to Fig. [1] it is in- 
teresting to note that the trend of increasing density near the 
core (R <10 kpc) based on the observed data appears to be 
one that crosses the two hydrostatic-equilibrium solutions, 
namely the LCD and HCD profiles. This suggests that, at 
least near the core of A2 199, deviations from such equilibria 
may be at play. In fact, I Sanders fc Fabianl (|2006T l hypothe- 
size the presence of a relatively weak (M ~1.5) isothermal 
shock in this region. The MHD simulations in 2-D planar 
geometry we performed here were intended to serve as a 
preliminary series of numerical experiments to assess pos- 
sible deviations from hydrostatic equilibrium in connection 
with transient AGN activity and/or remnant dynamical ac- 
tivity from the early formation of the galaxy cluster. To 
span a wide range of magnetic field effects we employed ini- 
tial strengths and topologies similar to those im posed by 
iParrish et all (120091 ) and lBogdanovic et al.l (|2009l '). The 2-D 
results now also form the basis for near-future 3-D simula- 
tions with MACH3. 

We were particularly interested in the initial response of 
the plasma near the core (t <1.3 Gyr, R =20 kpc) because 
radiative cooling and the HBI have been shown to lead to 
a runa way catastrophe in A2199 in this region in l ess than 
3 Gyr |Parrish et al.ll2009l ; IRuszkowski fc Ohll2010h . More- 
over, it has been proposed that mass and momentum trans- 



port from low-speed winds (~ 10,000 km/s) associated with 
AGN outbursts, can have a dramatic effect on the growth 
of the central supermassive black hole (|Ostriker et alil2010r i 
and, in turn, on the AGN-ICM coupling. We deliberately 
imposed in the initial thermodynamic state of the ICM a 
deviation from hydrostatic equilibrium, thereby producing 
subsonic expansion flow upon initiation of the simulation. 
By inducing a non-zero velocity field early in the simula- 
tion we allowed for (thermal) pressure work to contribute 
to the energy balance of the system. In fact, under the in- 
duced conditions, pressure work eventually overwhelmed the 
deleterious effects of the magnetic field, which in the ideal- 
ized hydrostatic case would have shut off the main supplier 
of heat to the region (in the absence of an AGN sources), 
namely thermal conduction. It is noted that only the scalar 
(thermal) pressure was included in this work; its importance 
in the presence of a subsonic velocity field and a magnetic 
field in this region suggests that the full pressure tensor m ay 
have to be accounted for as proposed bv lKunz et al.l (|2010T ll . 
In all three cases studied with different initial magnetic field, 
and for both initial density profiles (HCD and LCD), we ob- 
tained similar results by 2-D MHD simulations about the 
sustainment of the temperature near the core. 



5 CONCLUSIONS 

In this paper we presented results from a series of numerical 
experiments in two dimensions for the A2199 galaxy cluster. 
The simulations employed a relatively wide combination of 
initial conditions and simulation physics. It is acknowledged 
that the 2-D physical domains did not permit us to address 
a wide range of possible physics in clusters that are inher- 
ently 3-D, and the propositions presented herein shall be 
re-evaluated in three dimensions. Nevertheless, the number 
of cases performed in this study yielded constructive insight 
that now also forms the basis of our near-future simulations 
with MACH3. 

We chose A2199 largely because observed data for this 
cluster exist at a relatively small radius from its center 
(R ~1 kpc). As a consequence, our idealized 2-D axisym- 
metric simulations aimed to quantify the sensitivity of the 
ICM thermal balance on the plasma density profiles with 
less ambiguity compared to other clusters for which such 
data does not exist. The computed differences in the heat 
flux factor for the two profiles we used were not insignif- 
icant. Therefore, it appears it would be instructive that 
MHD parametric studies, performed to asses the sensitiv- 
ity of proposed dominant physics in the ICM and its sus- 
tainment/collapse (e.g. HBI, MTI, AGN heating, turbulence 
stirring), extended their range to include higher values of the 
plasma density near the cluster center. 

The results of the 2-D planar simulations demonstrate 
that imbalances in the system's hydrostatic equilibrium, in- 
ducing relatively weak flow dynamics, may overcome the ef- 
fects of the magnetic field on the diffusion of heat from the 
outer regions of the cluster. This marginalizes the signifi- 
cance of anisotropic thermal conduction over transient flow 
events in the ICM. Such events may possibly be associated 
with AGN pulsations at the cluster center, with characteris- 
tic dynamical times that are hundreds of times shorter than 
the Hubble time. 
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APPENDIX A: THE FULL MHD EQUATIONS 
OF MACH 

Al Conservation equations 

The MACH codes solve equations (|A"l"T) - (|A"6)) , the time- 
dependent, compressible, single-fluid, multi-temperature re- 
sistive MHD conservation laws, on a mesh that is composed 
of arbitrary hexahedral cells. 

dp 
dt 



dv 



= -V ■ (pv) 

= -pv ■ V« - V (p + Q + e R /3) 
+ J x B + pg 



(Al) 



(A2) 



p— = -pv ■ Ve e - p e V ■ V + V • («e ' VT e ) - <& eR 

+ J {rj-J)-J pCv {T e -T % )/r el (A3) 

en e 



p^ = -pv-Vei-(pi + Q)V -v + V -(Ri-VTi) 



OCR 
dt 



+ pCv (T e - Ti) /t. 



= —pv ■ X7e R — 4/3eflV ■ v 

+ V • (pXrVefl) + $ efl 
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^- = V x (v x B) - V x (fj. J) 

Vxl^Uvxf^ 

en e I V en e 
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In equations (|A1|) - (|A6|) p is the mass density of the fluid, 
v is the fluid velocity, B is the magnetic induction, t e n is 
the electron/ion specific internal energy, €r is the radiation 
energy density, and is related to the radiation temperature 
Tr via Stefan's constant a: €r = oTr. Stefan's constant is 
related to the Stefan-Boltzmann constant a = ac/4 where 
c is the speed of light in vacuum. T e u is the electron/ion 
temperature, p e /j is the electron/ion thermal pressure, Q 
is the artificial numerical compressional viscosity, J is the 
current density, n e is the electron number density, Cv is the 
specific heat, r e ; is the electron-ion equilibration time, e is 
the magnitude of the charge of an electron and fio is the 
permeability of free space. The electrical resistivity tensor is 
fj and R e /i is the electron/ion thermal conductivity tensor. 

Energy is conserved separately for electrons and ions. 
Models for radiation emission (thin limit), equilibrium 
conduction (thick limit) as well as non-equilibrium (flux- 
limited) radiation conduction also exist. The sesame tables 
serve as tabular equations of state for more than one hun- 
dred materials. The library is created and maintained by two 
groups at the Los Alamos National Laboratory (EOS and 
conductivity tables are the work of T-l; opacity tables are 
created and maintained by T-4). Although tabular values for 
the electrical and thermal conductivities are often available, 
the default models of our MACH version use the Spitzer- 
Braginskii conductivities parallel and perpendicular to the 
magnetic field separately for electrons and ions. Most of the 
physical processes are time-advanced in an implicit fashion 
that allows the Courant condition to be exceeded without 
exciting numerical instabilities. The spatial derivatives are 
obtained by integrating over finite volumes. The clean-up 
procedure applied to B to maintain its divergence at an ac- 
ceptably small level (solenoidal condition) uses a successive 
over-relaxation solver to iterate on the in-plane components 
toward V-i? = by solving for the scalar potential <f> (elliptic 
correction since a Poisson's equation is solved), and is b ased 
on the method proposed bv lBrackbill fe Barnes! (|l980t h 

In the last three decades MACH2 has been em- 
ployed to study a variety of plasma p roblems some of 
which include plas ma opening switches (|Buff et alj 1 19871 ; 



iDegnan et alj 1987f). gas and solid density z-pinch implo- 
sion physics l|Degnan et al.ll 19951 ; fShumlak et alj|l995h. ver y 
high-power plasm a source designs jMikellides et aljbooij) . 
magnetic nozzles (Mikellidcs et alJ 12002ft. and a variety of 
plasma thrusters lUMikellides et al.H2000l ; iTurchi et al.ll2000l; 
iMikellides fc Villarrealll2007l) . 

MACH is a three-temperature code, that is, the tem- 
perature of the radiation field Tr may be different from 
T e /i. In this case the radiation conduction equation (| A5|) 
is solved where both Planck (optically-thin limit) \ p and 
Rosseland (optically-thick limit) \r mean opacities must be 
supplied. The radiation coupling term is in general Q e R = 
acpx P (Te — Tr). Figure IAT1 compares the radiation emission 
rates ba sed on the MACH2 tabula r opacities and the for- 
mula bv lRvbicki fc Lightmanl (|l979h with a Gaunt factor of 
1.5 (note the factor ranges 1.1-1.5). 
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Figure Al. Comparison of the radiation emission rates 
based on the MACH2 ta bular opacities and the formula by 
iRvbicki fe Lightmanl jl979f ) with a Gaunt factor of 1.5. 



A2 Anisotropic Thermal Conduction Flux in 
MACH2 

In general, the conductive heat flux for electrons in a mag- 
netic field (in the absence the thermoelectric effect) consists 
of parallel, perpendicular and transverse components as fol- 
lows: 



1e = 1e// + 1e± + *? e A 

where 



(A7) 



1e// 
QeA 

and V 



= — «e//V//T e 
= — K e ±V±T e 

= -K fA ea x VT e 



// 



e B e B ■ V, Vj 



(A8) 

I — eses ) • V. The conduc- 



tivity coefficients are given by 



3.16re e o 

4.660,1 + 11.92 
Of + 14.79fii + 3.77 
fi e (2.5fie + 21.67) 



ReA - fi e + 14.7901+ 3.77 Ke0 ' (A9) 

and are the approximations o f lBraginskiil l| 19651 ) for the case 
of a fully-ionized plasma, with n e o — kBPeT e /m e . Thermal 
diffusion is solved for implicitly in MACH with an iterative 
Jacobi method. The rate of convergence of the Jac obi iter- 
ation may be increased by multigrid methods (e.g. iBrandtJ 
Il973h . which we have employed in the present MACH2 sim- 
ulations. 



